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Abstract 

At functional scales, cortical behavior results from the complex interplay of a large number of excitable 
cells operating in noisy environments. Such systems resist to mathematical analysis, and computational 
neurosciences have largely relied on heuristic partial (and partially justified) macroscopic models, which 
successfully reproduced a number of relevant phenomena. The relationship between these macroscopic 
models and the spiking noisy dynamics of the underlying cells has since then been a great endeavor. 
Based on recent mean-field reductions for such spiking neurons, we present here a principled reduction 
of large biologically plausible neuronal networks to firing-rate models, providing a rigorous relationship 
between the macroscopic activity of populations of spiking neurons and popular macroscopic models, 
under a few assumptions (mainly linearity of the synapses). The reduced model we derive consists of 
simple, low-dimensional ordinary differential equations with parameters and nonlinearities derived from 
the underlying properties of the cells, and in particular the noise level. These simple reduced models are 
shown to reproduce accurately the dynamics of large networks in numerical simulations. Appropriate 
parameters and functions are made available online for different models of neurons: McKean, Fitzhugh- 
Nagumo and Hodgkin-Huxley models. 

Author Summary 

Biologically plausible large-scale neuronal networks are extraordinary complex: they involve an extremely 
large numbers of neurons (on the order of hundreds of thousands), each of these having a complex, 
excitable nonlinear dynamics and subject to noise. This complexity motivated the use of deterministic 
low-dimensional descriptions of populations of neurons based on heuristic descriptions of the cortical 
function. These reduced models were successful in accounting for a number of phenomena in the brain 
such as spatio-temporal pattern formation. The relationship between well-accepted neuronal networks 
and heuristic macroscopic descriptions has since then been an important endeavor in the computational 
neuroscience community. We address here the problem of deriving such simple macroscopic descriptions 
rigorously based on the dynamics of large-scale neural networks. We apply our reduction method to the 
most prominent neuron models, and show that indeed, reduced models precisely account for macroscopic 
networks activity. By doing so, we can precisely compute the parameters involved in macroscopic models, 
and these are made available online for those neuroscientists aiming at precisely simulating large-scale 
networks of McKean, Fitzhugh-Nagumo or Hodgkin-Huxley neurons. Analytical developments are also 
provided in a simplified model to illustrate the reduction method. 
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Introduction 

The activity of the brain is characterized by large-scale macroscopic states resulting of the structured 
interaction of a very large number of neurons. These macroscopic states correspond to the signal experi- 
mentally measured through usual recording techniques such as extracellular electrodes, optical imaging, 
electro- or magneto- encephalography and magnetic resonance imaging. All these experimental imaging 
protocols indeed record the activity of large scale neuronal areas involving of the order of thousands 
to millions of cells. At the cellular level, neurons composing these columns manifest highly complex, 
excitable behaviors characterized by the intense presence of noise. Several relevant brain states and func- 
tions rely on the coordinated behaviors of large neural assemblies, and resulting collective phenomena 
recently raised the interest of physiologists and computational neuroscientists, among which we shall cite 
the rapid complex answers to specific stimuli [I], decorrelated activity (2[[3], large scale oscillations ^, 
synchronization |5], and spatio-temporal pattern formation (6|[7j. 

All these evidences motivate the development of models of the collective dynamics of neuronal pop- 
ulations, that are simple enough to be mathematically analyzed or efficiently simulated. A particularly 
important problem would be to derive tractable macroscopic limits of the widely accepted and accurate 
Hodgkin-Huxley model [sj. However, describing the activity of a network at the cellular scale yields 
extremely complex, very high dimensional equations that are mathematically intractable and lead to 
excessively complex and time consuming numerical simulations. 

The question of the macroscopic modeling of cortical activity and their relationship with microscopic 
(cellular) behavior has been the subject of extensive work. Most studies rely on heuristic models (or 
firing-rate models) since the seminal works of Wilson, Cowan and Amari [9 10 . These models describe 



a macroscopic variable, the population-averaged firing-rate, through deterministic integro-differential 
or ordinary differential equations. Analytical and numerical explorations characterized successfully a 



number of phenomena, among which spatio-temporal pattern formation and visual illusions (see 11 for 
a recent review). This approach was complemented by a number of computational studies introducing 
noise at the level of microscopic equations, the effect of which vanishes in the limit where the number of 
neurons tends to infinity. These approaches are generally based on simplified neuron models and make 
strong assumptions on the dynamics (e.g. sparse connectivity |12| , Markovian modeling of the firing 
and van Kampen expansion [l3]). Relationship between spiking neuronal networks and mean firing rates 
in simplified models and deterministic settings has also been the subject of a number of outstanding 



researches 14 15 . These averaging techniques were based on temporal averaging of periodic spiking 



behaviors. For instance, in 15 , the author presents a reduction to Wilson-Cowan systems for the single- 
cell deterministic Morris-Lecar system, taking advantage of the separation of timescales between slow 
synapses and cell dynamics. In contrast with these researches, we propose a mixed population and 
temporal averaging for stochastic networks, taking advantage of the collective effects arising in large 
networks. 

Despite these efforts, deriving the equations of macroscopic behaviors of large neuronal networks from 
relevant descriptions of the dynamics of noisy neuronal networks remains today one of the main challenges 



in computational neurosciences, as discussed in P. Bressloff's review 11 . This is precisely the aim of 



the present manuscript. Inspired by statistical mechanics methods, we start from rigorously derived 



limits of neuronal network equations 16 with excitable dynamics of Hodgkin-Huxley type. The resulting 
equations, referred to as the mean-field equations, are hard to interpret and to relate to physical observable 
quantities. In the gas dynamics domain, mean-field equations such as Boltzmann's equation were used 
to derive the behavior of macroscopic quantities such as the local density, macroscopic local velocity and 
local temperature fields, in relationship with the microscopic activity of the particles, and resulted in 
the derivation of the celebrated Navier-Stokes equations that provide important information on the fluid 
dynamics. In our biological case, a particularly important quantity accessible through measurement is 
the macroscopic variable corresponding to an averaged value, over neurons at certain spatial locations, 
and on a specific time interval, of the activity of each cell. By doing so, we will reduce the complex high 
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dimensional noisy dynamics of microscopic descriptions of neural networks into a simple, deterministic 
equation on macroscopic observables. 

In the Material and Methods section, we will introduce the basic network equations and their mean- 
field limits, and describe the methodology we propose for deriving macroscopic equations. This method 
reduces the dynamics of the average firing-rate to the knowledge of a particular function, the effective 
non-linearity, which can be numerically computed in all cases. This methodology is put in good use in the 
Results section in the case of the McKean, Fitzhugh-Nagumo and Hodgkin-Huxley neurons. In each case, 
the effective-nonlinearity is numerically computed for different noise levels. The reduced low-dimensional 
macroscopic system is then confronted to simulations of large networks. The discussion section explores 
the limitations and some implications of the present approach. In the appendix, analytical developments 
are also provided in the case of the deterministic McKean neuron in order to illustrate the mechanics of 
the reduction. 



Material and Methods 

In this section we introduce the original networks models, their mean-field limits and the formal derivation 
of the dynamics of averaged firing-rate models. This approach will be used in the result section to derive 
macroscopic limits and demonstrate the validity of the reduction. 

Neurons and networks 

We consider P populations networks composed of {Na)a£{i---p} neurons. Functional structures involve 
very large numbers of neurons Na in each populations, and similarly to the case of gas dynamics, we will be 
interested in the limit where all Na tend to infinity (the mean-field limit) in order to exhibit regularization 
and averaging effects. Each neuron i in population a is described by the membrane potential vl (also 
called activity of a neuron in this paper) and additional variables gathered in a d-dimensional variable 
Zl , representing for instance ionic concentrations in the Hodgkin-Huxley model, or a recovery variable in 
the Fitzhugh-Nagumo or McKean models. These variables satisfy a stochastic differential equation: 



dvl = (Paiv', Z') + I"{t) + Jijvl *h]dt + aadWl 
dZl ^Ga{vl,Zl)dt + Ta dBl 



(1) 



In this equation, the functions Fa and Ga describe the intrinsic dynamics of all the neurons in population 
a. The parameters CTq, G M and Fq, S W^'^^^ describe the standard deviation of the noisy input received 
from each component of X"^ — (w*,Z'), classically driven by independent one (resp. d) dimensional 
Gaussian white noise {WD (resp. El). I°'{t) represents the external input neurons of population a 
receive. The coefficients Jy are the synaptic weights of the connection j — >■ i. Spiking interactions 
between neurons are here modeled as the convolution of the presynaptic membrane potential with the 
impulse response of the synapse noted /i, where h{t) = -^e^^Ltyo with Tg = 10 ms is the characteristic 
time of the synapses and lt>o the Heaviside function. This is an idealization synapses dynamics under 
the assumption that the synapse is linear. 

Three main models used in computational neuroscience addressed in the present manuscript are the 



McKean model 17 , the Fitzhugh-Nagumo model \18^ and the Hodgkin-Huxley model These models 



are parametrized so that the time unit is in milliseconds. 



The McKean model is a piecewise continuous approximation of the Fitzhugh-Nagumo model allowing 
explicit calculations ^17] . The membrane potential of neuron i is written and its adaptation variable 
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w*. They satisfy the following network equations: 



dvl 



with 



dwl — ^wivl — wl + b)dt 

-Iv — (I + c)a if t; < —a 

/(u) — cv ii a < V < a 

-Iv + {I + c)a if a < w 



(2) 



and e„, ~ 0.1, / = 1, a = 1, c = 0.5, and h — 0. 



The Fitzhugh-Nagumo model was introduced in 18 and has been widely studied as a paradigmatic 



excitable systems for its simplicity and ability to produce spiking behaviors. This model was introduced 
in order to allow analytical study in a simpler model mimicking the essential features of the Hodgkin- 
Huxley neuron. Fitzhugh-Nagumo models describe the activity of the membrane potential and recovery 
variable of neuron i in the network through the equations: 



dvl-- 
dwi 
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I' + J'^j i^i *h)]dt + a'dW, 



b)dt 



(3) 



with e„ = 0.08, a = 0.8 and b = 0.7. 



The Hodgkin-Huxley model: Probably the most relevant neuron model from the electrophys- 
iological viewpoint, the Hodgkin-Huxley model describes the evolution of the membrane potential in 
relationship with the dynamics of ionic currents flowing across the cellular membrane of the neuron, de- 
veloped in [s] after thorough observation of the giant squid axon revealed the prominent role of potassium 
and sodium channels for excitability, and leak chloride currents. Networks of Hodgkin-Huxley neurons 
are described by the stochastic differential equation: 



Cdvi ^ 



dnl 



where 



a„(w) = 0.01- 



dm\ 
dhi 
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- gKin')Hvl - Ek) - gNa{m'fh\v, 

-gL{vl~-EL) + Y.,'h {vl*h)yt 
'an{vl){l~n')-pn{vlW)dt 
'am{vl){l - m*) - I3m{vl)m')dt 
\ah{vl){l ~ h') ~ Ph{v\W)dt 



t 

a'dW, 



En a) 

'i 



exp(^) ~ 1 



a-miv) = 0.1- 



25 



exp(^) - 1 



0.125 exp( — ). 



4.exp( — ) 



a^(z;) = 0.07exp( — ) 



exp( 



10 ' 



1 



(4) 



ancJ^C = IfiF/cm'^, gx = 36mS/cm'^, g^a — 12QmS/cm?, gL — 0.3mS'/cm^, Ek — —12mV, E^a = 
120my and El = lO.GmV . The dynamics of this system shows deep non-linear intricacies even in the 
case of deterministic, single-neuron system. 



•^Parameters are adjusted so that the resting state of the membrane potential v is equal to OmV. 
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Mean-field limits 

These network equations are extremely complex to analyze and simulate. However, in the mean-field 
limit, one can access to the asymptotic behavior of neurons in large networks. We consider that synaptic 
weights are heterogeneous and have a typical mean and standard deviation depending on the populations 
they belong to, e.g. Jij ^ J\f{^^, -^) where Ja-, is the typical connectivity weight between populations 
a and 7 and A the synaptic weights disorder. In this case, the network is described by a set of delayed 
stochastic differential equations, and the propagation of chaos applies, as shown in [16] . This means that 
all the neurons in the same population (which we call a) follow the same law. This law is described by 
the following mean-field equation 



Observe that this equation, describing the process governing each neuron in population a, is very similar 
to the network equation ([T]). Only the interaction term is different since it involves the computation of 
an expectation of the activity. The mathematical study of this type of equation is generally extremely 
complex: in particular, it is not a stochastic differential equation, but rather an implicit equation in 
the set of stochastic processes. In our present approach, we will manage to bypass this difficulty by 
considering the macroscopic activity of a population. 

Firing Rates, Macroscopic Activity and Dynamics 

The probabilistic behavior of neurons in the limit N ^ 00 being characterized, we can define macroscopic 
variables caracterizating the state of the population. 

The firing-rate is usually considered a good candidate as a macroscopic descriptor of the population 
activity. Heuristically, its corresponds to the average number of spikes fired in a certain time window. 
However, counting spikes is a non-linear and non trivial operation which is difficult to handle mathemat- 
ically. 

We will rather use a simpler definition of the macroscopic activity: the average within a population 
and during a certain time of the membrane potential of neurons. This definition has the interesting 
property of being a linear transformation of the activity of neurons. It is also closely related to the firing- 
rate. Indeed, given that neurons communicate via spikes which are stereotyped electrical impulses of 
extreme amplitude, integrating the value of the membrane potential during a time window and dividing 
by the area under a spike provide a rough estimate of the number of spikes emitted. However, the two 
measures are not exactly similar: this definition of macroscopic activity is impacted by the activity under 
the firing threshold of neurons as opposed to pure firing-rates. 

More precisely, we define the macroscopic activity of population a as the spatial mean of the 
membrane potential over the neurons in the populations and temporally convolved with the "time 
window" function g of width = 100 larger than the duration of the spikes, but small enough to 
resolve fine temporal structure of the network activity. The propagation of chaos property ensures that 
all neurons are independent; this allows to identify the population-averaged voltage with the statistical 
expectation of the voltage variable in the mean-field limit and leads to the following definition of the 
macroscopic activity: 




a 



(5) 




(6) 



•^Concretely, g{t) 



1^ 
K 



_ t~ 

e with s chosen so that Kg(d/2) = 0.01 and K so that J^ g{t)dt = 1. 
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Thanks to the hnearity of synapses, using this definition together with the mean field equation ([5]) leads 
to the following caracterization of the macroscopic activity. 

p 

= E *g + Y, Jo.p{y^ * h) + /"(t) (7) 

,3=1 

where X" = (w", Z"), /" denotes /" * g. 

This equation is not closed because of the term E * g which is not expressed as a function 

of the macroscopic activity. Because of the nonlinearity of F^, it is not likely that this quantity only 
depends on i/". Moreover, the membrane potential depends on additional variables and the macroscopic 
activity hence involve expected value of functions of these variables convolved with the time kernel g. 

Computing the effective non-linearity 

We now attempt to make sense of the term E [Fa{Xf )\ * g to find a closed formula on the macroscopic 
activity. We introduce the ansatz that this term can be written as the sum of a linear functional of the 
macroscopic activity and a non-linear term applied to the effective input to a population (including the 
synaptic connections). Defining this effective input x" as: 

p 

x"'!/^ J,^(j.''*/i)+/"(0 (8) 

/3=1 

the ansatz reads: 

E[F„(Xf)]*.g = i(i.") + 5(^") (9) 

where L is a linear functional and 5 is a non-linear mapping which remains to be determined. The 
choice of this ansatz was motivated by the analytical treatment of networks of deterministic McKean 
neurons (see appendix), which naturally exhibits a relation like ([£]). The choice of the linear functional 
L is dictated by the neuron model used. In general, the function S has to be computed numerically. 

To evaluate S, one need to assume that both the inputs and the synapses are slow compared to the 
dynamics of the neurons. Thus, the effective input can be considered as constant during the time 
the neurons reach an asymptotic regime related to that input state. If this regime is stationary, the 
value function S at will be the average of i^^ — L applied to that stationary stochastic process, and 
therefore will provide a quantity only depending on x". If the regime is periodic in law, then taking 
a time window g of size 9 larger than the periotQ will also yield a constant value for our nonlinear 
term. Because the mean-field equation ([5| and a single neuron equation only differ in the interaction 
term which is assumed constant here, the computation of S{x°') simply corresponds to computing the 
temporal average of Fa — L along the trajectory a single noisy neuron forced with a constant input 
. Actually, for readability reason^ we will rather focus on the function S which is simply defined as 
Six"') = S{x°')+x°'. 

One of the pitfalls of this methods occurs if the behavior of the mean-field equation depends on the 
initial condition and present multiple stable stationary or periodic attractors, the quantity E [i^Q,(X")] * g 
can take different values depending on the initial condition. A neuron model (together with a particular set 
of parameters) will be said to be of regime I if there is only one of these attractor for any initial condition. 
Similarly, if there are p-attractors, the neuron model is said to be of regime p. Figure [T] shows different 
values for S{x°') when starting from different initial condition. Figure Ilia) shows the solution for regime 
I and figure [l]b) for regime II. The points obtained (see Fig. [T]) are then segmented into a few clusters 

■^we recall that 8 was assumed to be small enough to resolve the temporal structure of the neuron's activity. 
^This is motivated by the combination of equations (ml, JSjl and (j9l into equation ||10[l 
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(in our case, one or two in regimes I and II respectively) and smoothed out into a surface (or a union of 
surfaces) see figure |3] This procedure is relatively time consuming. The result of extensive simulations on 
the McKean, Fitzhugh-Nagumo and Hodgkin-Huxley models, using this numerical procedure, are freely 
available online as well as the algorithm generating these data. 




(a) Regime I neuron (b) Regime II neuron 

Figure 1. Value of S{x°') computed for 10 different initial conditions (and 200 in the inset of (b)) for 
the Hodgkin-Huxley model with noise cr = 2 (a) and a — (b). This shows that the level of intrinsic 
noise change the regime of a neuron. 



Linear part L for the different models 

Identifying L for a given model, consists in gathering the most linear terms in the intrinsic dynamics of a 
single neuron F. This linear terms can be time-delayed. A careful study of the different models suggests 
that the linear parts for the different neurons can be chosen as in table [T] 



Model 


Linear part 


McKean 




Fitzhugh-Nagumo 




Hodgkin-Huxley 





Table 1. Linear part L for the different models. 5 is the Dirac function centered at and 
hr^. = eu)e~'^™*lt>o for the first two models. 

For McKean and Fitzhugh-Nagumo models, the convolution accounts entirely for the existence of the 
adaptation variable w which will therefore still be present in the reduced model. The choice of the linear 
part of the McKean model can be understood better by reading the appendix where it naturally comes 
out of the computations. The choice of the coefficient 4/3 for the Fitzhugh-Nagumo model relies on an 
analogy to the McKean neuron. It corresponds to the (absolute value) of the slope of the straight line 
approximating the negative decreasing part of the non-linearity Fitzhugh-Nagumo v — As for the 
Hodgkin-Huxley, we were not able to extract a linear delayed term in the intrinsic dynamics of a neuron 
so we simply chose the linear decay already present in the equations. 
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Simulation of the macroscopic equations 

For regime I neurons, when the effective non-hnearity S is univalued, simulations of the macroscopic 
equations simply reduce to solving numerically the following ordinary differential equation 

p 

= L{iy") + 5( ^ J^p{,^^ * h) + /"(O) (10) 

13=1 

For regime II neurons, when the initial condition is not in the bistable region, we will consider that the 
averaged system pursues on the initial attractor (fixed point or spiking cycle) when possible, and switches 
attractors if the activity brings the system in regions where the initial attractor disappears. In details, 
let us denote by S{x°', 1) the branch of stable fixed points, defined as long as x" < Ih, and by S{x°',2) 
the branch of periodic orbits defined for x" > Iflc- The macroscopic activity of population a in a 
P-population network hence satisfies the equations: 

= L{iy^) + Jo.pi'^^ * h) + /",p(t)) ^^^^ 

[p{t) = 4»-/„(5p,i - 6x''-ipLc^p,2- 



Results 

In this section, we evaluate the accuracy of the reduced model presented above for the three neuron 
models McKean, Fitzhugh-Nagumo and Hodgkin-Huxley. First, we address the computation of the 
effective non-linearity both in the deterministic and noisy cases. Second, we confront the time course 
of the macroscopic activity calculated according to our reduction against the a posteriori average of the 
activity of a spiking network. 



Effective non-linearity without noise 

In the deterministic reduction (small noise limit), the effective non-linearity can be understood through 
the bifurcation analysis of a single neuron. Indeed, the effective sigmoid amounts to computing the 
temporal average v of the voltage solution of the deterministic single-cell system upon variation of the 
input. The result of that analysis using the numerical software XPPAut 19 is displayed in Figure Fig. [2] 
In these diagrams, we colored regions of stationary solutions (green), periodic solutions (purple) and 
bistable regimes with co-existence of a stable stationary solution and a stable periodic orbit. 



The McKean neuron 

Although the deterministic McKean neuron is analytically treated in the appendix, we now consider 
it under the angle of bifurcations for consistency with the other models. In the McKean neuron, the 
non-differentiable, piecewise-continuous nature of the flow gives rise to a non-smooth Hopf bifurcation 
associated with a branch of stable limit cycles. The emergence of the cycle arises through a non-smooth 
homoclinic bifurcation, hence corresponding to the existence of arbitrarily slow periodic orbits, typical 
of a class I excitability in the Hodgkin classificatior{^ In this model, an important distinction is the 
absence of bistable regime: the average variable has a unique value whatever the initial condition (except 
singularly chosen equal to the unstable fixed point) and whatever the input chosen. In the present case, 
stable permanent regimes are unique. Therefore, there exists a single-valued function S defining the 



dynamics of the macroscopic activity through equation ( 10 ) 



"Observe that Hodgkin classification is different from the classification in regime we propose in this paper 



9 






H 














FLC 


FLC 














H 





PLC 










* H 



5 10 15 20 

I 



I I 

(a) McKean bifurcation diagram (b) FitzHugh-Nagumo bifurcation (c) Hodgkin-Huxley bifurcation dia- 

diagram gram 



Frequency 




Frequency 




Frequency 

0.1 





1 


1 ^^^..---'-'''^ 






H ^ 


FLC ^ 













(d) McKean frequencies 



(e) FitzHugh-Nagumo frequencies 



10 I 20 

(f) Hodgkin-Huxley frequencies 



Figure 2. Bifurcation diagrams (upper row) representing the temporal average of the solutions (i.e. 
the fixed points and average value in the case of periodic orbits) and frequency of the regular spiking 
regime (lower row) in the McKean model (left), Fitzhugh-Nagumo model (center) and Hodgkin-Huxley 
model (right). 
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The Fitzhugh-Nagumo model 

The bifurcation diagram of the Fitzhugh-Nagumo neuron as a function of the input level (see figure [2(b)| ) 
presents multi-stability. For small negative input, the system presents a stable fixed point. Increasing the 
value of the input makes the fixed point lose stability through subcritical Hopf bifurcation, and unstable 
limit cycles appear. These limit cycles undergo a fold and a branch of limit cycles appear, overlapping 
in a small parameter region the state where a stable equilibrium exists. This small parameter region 
again corresponds to a bistable regime with co-existence of a resting state and of a regular spiking 
behavior. This branch of stable limit cycles corresponds to a regular spiking regime. As the input is 
further increased (in a biologically unplausible range), the same scenario arises symmetrically: the branch 
of stable periodic orbits undergoes a fold of limit cycles, a branch of unstable periodic orbits emerges 
from this bifurcation and connects with the unstable fixed point at a subcritical Hopf bifurcation and 
the unstable fixed point gains stability. Here again, the neuron corresponds to a class II excitability and 
a regime II, but the system could be well approximated by a excitability 1/ regime I since the bistability 
region is of very small extent and the periodic orbits have very small periods when appearing, close from 
a class I excitability. Thus, we consider that the macroscopic activity evolves according to equation ( 10 1 
as illustrated in the following. 



The Hodgkin-Huxley model 

This model displays qualitatively the same behavior as the Fitzhugh-Nagumo model but with significant 



quantitative differences: the bifurcation diagram of the Hodgkin-Huxley neuron (see figure 2(c) ) also 
displays multi-stability. We observe the existence of a branch of stable fixed points (green region) that 
undergo subcritical Hopf bifurcation for = Ih, associated to a family of unstable limit cycles (pink 
dotted lines represent the average value of v along the cycle) . This family of unstable limit cycles connects 
with a branch of stable limit cycles (pink solid line) through a fold of limit cycles bifurcation. These 
stable limit cycles are the unique attractor for large input (implausibly large input values will nevertheless 
see these cycles disappear in favor of a high voltage fixed point). The system presents a bistable regime 
(yellow input region) where a stable fixed point and a stable periodic orbit co-exist. The frequency along 
the cycle (Fig. |2(f)[ ) shows a class II excitability in the Hodgkin classification: oscillations appear with a 
finite period and a non-zero frequency. The hysteresis present in the yellow region corresponds to what 
we called a regime II. In simulations, when the initial condition is not in the bistable region, we will 
consider that the system pursues on the initial attractor (fixed point or spiking cycle) when possible, and 
switches attractors if the activity brings the system in regions where the initial attractor disappears, as 
explained in the Material and Methods section. This method is chosen here because the bistability only 
appears for small noise, regimes in which switches between the different attractors are rare. 



Effective non-linearity with noise 

When considering noisy input, no analytical approach could be performed. We use the method described 
in Material and Methods and plot the results in Figure [3] where the effective non-linearity is displayed as 
a surface obtained as a function of noise and effective input. It appears relatively clear in the figure that 
noise tends to have a smoothing effect on the sigmoids. This effect is particularly clear in the Hodgkin- 
Huxley model where a multivalued function (regime II) is turned into a single valued smooth function 
(regime I). In the case of the Fitzhugh-Nagumo model, we observe that the regime II is not observed in 
simulations in the presence of noise. This is due to the smallness of the parameter region corresponding 
to the bistable regime, and the averaged system can be well approximated by regime I dynamics. In the 
case of the Hodgkin-Huxley network, there are clearly two different behaviors depending of the level of 
intrinsic noise. When the noise is small (resp. large) the neuron is regime II (resp. I). Interestingly, this 
shows how a strong noise can qualitatively simplify the macroscopic dynamics of a network. 
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a) McKean b) Fitzhugh-Nagumo c) Hodgkin-Huxley 




Figure 3. Effective non-linearities in the McKean, Fitzhugh-Nagumo and Hodgkin-Huxley model. 
Observe that noise tends to have a smoothing effect on the sigmoids. For the Hodgkin-Huxley model, we 
have empirically chosen a noise threshold under which the neuron was considered regime II and above 
which it is regime I. There are thus 2 branches below the threshold and only one above. 



Comparison between reduced model and averaged spiking network 

We now simulate large networks of McKean, Fitzhugh-Nagumo and Hodgkin-Huxley neurons, compute 
numerically their macroscopic averaged activity, and compare the dynamics of this variable to simulations 
of the reduced ordinary differential equation involving our effective non-linearity function in all three cases. 
We consider that there are P = 5 populations of = 200 neurons (i.e. — Np for all a, All the 



neurons in the same population receive the same input corresponding to one color in figure 4(a) Some 
parameters are constant for all simulations: — 10ms, 9 — lOOms. Simulation dependent parameters 
are detailed in the caption of figure [4] which gathers the comparison for the different models. 

The comparison for the McKean and Fitzhugh-Nagumo neurons show a precise match see figures 
|4(b)]|4(^ and [4(d)] even when the strength /i = 1 of the connections is strong enough to significantly 



modify the shape of the input signal. Note that, although the mathematical derivation assumed slow 
inputs and synapses, these simulations validate the reduction even with faster synapses (r^ = lOms) and 
inputs evolving at the scale of hundreds of milliseconds. Also observe that adding frozen noise on the 
connections with A = 1 in figure [4(d)] does not significantly change the accuracy of the match. 

The comparison for the Hodgkin-Huxley neuron are only a relative success. First, the reduced model 



is only accurate when the connections between populations are kept small, i.e. fx = 0.1 in figure 4(e) 
This is why the simulations for the Hodgkin-Huxley networks show patterns very similar to the inputs 
When the connection strength increases the quality of the match significantly decreases (not displayed 
here). Second, the algorithm proposed in Material and Methods to simulate the regime II networks fails 
reproducing faithfully the averaged spiking network, see figure [4(f)[ More precisely the jump from one 
branch to the other is not well approximated. We observe that the populations having fast switches in 
the multivalued region (blue and green) or that do not intersect the multivalued input region (cyan and 
purple) are precisely recovered. However, the red population, spending much time in the multivalued 
region, is not well approximated by the macroscopic activity model: the network equations may randomly 
switch from spike to rest, which produces the irregular macroscopic activity, whereas the smooth firing 
rate does not show such switches. Notwithstanding these unavoidable errors, we observe a fair fit of the 
macroscopic activity model, which recovers most of the qualitative properties of the network activity. 



Discussion 



Even if collective phenomena arising in large noisy spiking neural networks are extremely complex, we 
have shown that the macroscopic activity of the network can be consistently described by simple low 
dimensional deterministic differential equations. The parameters and non-linearity involved are deter- 
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(a) Inputs (b) McKean, regime I. cr = 0.1,^ = (c) Fitzhugh-Nagumo, regime I. cr = 

1,A = 0.5,Ai=l,A = 




time (ms) time (ms) time (ms) 

(d) Fitzhugh-Nagumo, regime I. cr = (c) Hodgkin-Huxley, regime I. tr = (f) Hodgkin-Huxley, regime II. u = 
0.5,Ai = l,A = l 1.5,At = 0.1, A = 0.1, = 0.1, A = 



Figure 4. Comparison between tiie direct computation of tlie averaged macroscopic variables 
computed through direct simulation of the network equations (plain lines) and simulations of the 
macroscopic equations (dashed lines) for networks of P = 5 populations with TVq, = 200 neurons 
(different populations are depicted in different colors). The connections are drawn randomly according 
to the Material and Methods description and the random seed is the same for simulations 4(c) and |4(d)] 
The inputs J to the McKean and Fitzhugh-Nagumo networks are shown in (a). The inputs to the 
Hodgkin-Huxley networks correspond to 100/ + 10 where / is shown in (a). The simulations where done 
using a stochastic Euler algorithm with T = 15000 (resp. T — 30000) timesteps of size c?i = 0.1 (resp. 
dt — 0.05) for McKean and FitzHugh-Nagtimo (resp. Hodgkin-Huxley) networks. 
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mined by the type of neurons eonsidered. Depending on the neuron model the non-hnearity can be a 
well-behaved function (which we call regime I) or a more complicated multivalued function (which we call 
regime II in case of two values). The three neuron models we considered (McKean, Fitzhugh-Nagumo 
and Hodgkin-Huxley) are regime I when the intrinsic noise is strong in the network. However, with weak 
noise the Hodgkin-Huxley model is regime II, in which case the low-dimensional model proposed is more 
complicated (involving jumps between attractors). Comparisons of the averaged dynamics of spiking 
networks with the reduced equations showed a very precise fit, even for initial conditions independent of 
the network' s initial conditions, for regime I neuron models. However, for the regime II neuron models, 
the reduction accuracy is not as good. Indeed, noise will induce random switches from one attractor to 
the other, which cannot be handled through reduced methods. Yet, the reduced model recovers the main 
qualitative features of the signal, but in the bistable regions, quantitative distinctions arise. 

The reduction accuracy is significantly better for McKean and Fitzhugh-Nagumo neuron models 
than for the Hodgkin-Huxley model. Indeed, the reduction for the latter becomes irrelevant for strong 
connections between neurons whereas it is not the case for the former. We believe this is not due an 
inherent difference between the models, but rather an incapacity from us to determine a relevant linear 
part L for the Hodgkin-Huxley model. Indeed, as shown in tablejlj there is some time-delayed information 
in the linear part of McKean and Fitzhugh-Nagumo whereas there is simply a linear instantaneous decay 
for Hodgkin-Huxley model. Further work on reducing Hodgkin-Huxley networks should include a way to 
add information in the linear part of this model. 

This reduction holds on a number of assumptions imposed by the mathematical approach: we need 
connections within populations to be homogeneous and scaled by the number of neurons in the population, 
which has to be large enough for averaging effects to occur (i.e. for the mean- field reduction to hold). 
Besides, the inputs to the neurons in the same population are considered identical. More importantly, 
the reduction is largely based on the linearity of synapses. Although the slowness of synapses and inputs 
were required for the mathematical derivation, simulations have shown that the reduction was quite 
robust to increased speed for synapse and inputs. It is also fair to observe, that this reduction can not 
be applied directly to real biological tissues. All the neurons are assumed identical and the propagation 
of action potential along axons (which is important in measurements) has been ignored. In the future, 
an extensive work involving detailed biological knowledge and heavy computer simulations could make 
possible to apply this reduction to biological neural assemblies. 

It is important to note the huge complexity reduction obtained: in the case of Hodgkin-Huxley 
networks, we reduced a system of AP N stochastic differential equations with P populations into a system 
of P deterministic, ordinary one dimensional differential equations. The nonlinear transforms computed, 
as well as code for the simulations, are provided freely online. For efficient simulation of large-scale 
neuronal spiking networks with noise, if one is interested in computing the mean macroscopic activity, 
simulating the reduced model is a precise and simple choice that shall be considered for efficiency. 



This study quantifies the stabilization properties of the noise, that were already discussed in 20 
controlling the shape of the effective non-linearity of the reduced model. Noise tends to act as a linearizer: 
when the noise is strong, the network macroscopic activity tend to evolve more linearly. The example of 
Hodgkin-Huxley model shows it can even change a neuron model from regime II to regime I. This implies 
that knowing the value of the intrinsic noise in biological tissues could be a good indicator to evaluate 
their level of non-linearity. 
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Appendix A: Analytical derivations for the deterministic McK- 
ean model 

Using a McKean network ^ in the mean field reduction ^ in the case cr* = leads to the reduction 



dv° 
w° 



_ + ^-^^ j^^E [v-r] * h + I'^it) ]dt 



(12) 



The imphcit integration of the adaptation equation in ( 12 ) shows that w" ~ v" * hr^ + b, where 
hT^{t) = et„e~'™lt>o. Using the conimutativity of the convohition and equation ([6]), we obtain the exact 
macroscopic equation for the McKean neuron: 



(13) 



The only unknown term in the formula above is E(/(w")) * g. To compute this term, consider system 



(12 1 where = J^-y^i Ja-y^ b7] * h + I^it) is assumed to be constant. Indeed, as mentioned before, 



the following derivation relies on the fact that the inputs and the synapses are so slow that the adiabatic 
assumption holds (i.e. reaches its equilibrium value immediately). Thus, can be considered 
uncoupled from the others populations and it has the same dynamics as a single McKean neuron, whose 
phase plane is shown in figure|5] Under the assumption that the recovery variable is very slow (e^, <^ 1), 




Figure 5. Phase plane of the deterministic McKean neuron (equations ( 12 ) with 

x°' = X)/3=i U"*^ + I^it) a constant). When the blue line and any of the decreasing green line intersect, 
then there is a stable fixed point. When the blue curve intersect the increasing green line, then there is 
a periodic orbit (this is the case shown here). The periodic orbit corresponds to the non-smooth orange 
trajectory composed of two branches on the slow manifold (single-arrowed segments) and horizontal 
double-arrowed segments correspond to the fast transitions. 



the state of neurons in population a is essentially projected on one of the two slow manifold, corresponding 
in the phase plane Figure [5] to the single-arrowed orange branches of the w-nullcline. Extemely fast 
switches between these two branches of the slow manifold occur when the trajectories reach an extremity 
of the manifold. Except during the very fast transitions, it holds that f{v°') — —Iv" ± {I + c)a + x" and 
hence 

E [/(u")] -m + E(±)(/ + c)a + 
with E(±) dv"{t) - dv°'{t) = V(v°' > 0) - P(w" < 0). Therefore we have: 

E(/(w")) * g = -Iv'^ + + c)a{V{v°' > 0) - P(i;" < 0)) * .g + 
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We write S'(x") S K the value toward which the function t ^ [^(P(i;" > 0) - P(u" < 0)) * gj{t) 

converges when the McKean is stimulated by a constant input x" (thus we discard the initial transient). 
Computing P(w" ^ 0) * g for a constant effective input x°' amounts to computing the proportion of time 
system a McKean neuron spends on (or close to) the slow manifolds = —Iv" ± (/ + c)a + x°'. This 
can be performed analytically. Indeed, there are two different cases: 

• If a;" < -(1 -c)a + b (resp. > (1 - c)a + b). 

Then the system has a single stable fixed point on the negative (resp. positive) slow manifold. In 
figure [5j this corresponds to the blue curve crossing the green piecewise cubic where the latter is 
decreasing. In this case, S{x°') — —1 (resp. S{x°') = 1). 

• If -(1 -c)a + b< < (1 - c)a + b. 

Then the system is oscillating on a deterministic limit cycle represented in orange in figure [5] In 
this case, S{x°') = j'+[x'')+t-[x'^) where T~{x") (resp. T+(a;")) is the duration it takes for the 
system to go along the negative (resp. positive) part of the slow manifold. Following we can 
access these values. Indeed, assume the fast membrane potential immediately goes to one of the 
slow nuUclines. This gives the equation: —Iv"" ± (Z + c)a — + a;" = 0. Injecting this in the slow 
equation and integrating along relaxation orbit (orange path in figure [5| leads to 

I r"^"° dw I , f{l + 2c/l + c)a + b-x° 



J-ca+x" + + [l + c)a + x"" + lb twil + l) V (l-c)a + &-x" 
Similarly, 

I /(l + 2c/Z + c)a + a;"-6 



^ ' e,.,{l + l) 1-ca 



e„(l + (l-c)a + a;^ 

Therefore, for e] - (1 - c)a + 6, (1 - c)a + b[ 



((l+2c/J+c)a+6-2;°) ((l-c)a+i:° -b) \ 
((l-c)a+6-rE = ) ((l+2c/(+c)a+a;°-6) / 

^(x") = ^ ^ (14) 

((l+2c//+c)a+b-2;°) ((l+2c/(+c)a+i:° -b) \ 
((l-c)a+b-a;") ((l-c)a+i;°-b) j 



log 
log 



This function is shown in figure |6] It is a non-smooth sigmoidal function with vertical tangents at 
— (1 — c)a + b and (1 — c)a + b. This corresponds to the transition from a fixed point to the oscillatory 
pattern. Note that it is identical to the bifurcation diagram of Fig. [2(a)] 



Based on equations ( |13[ ) and the definition of the sigmoid ( 14 ) , we are now in position to define an 
averaged model describing the evolution of the macroscopic population activity. It takes the form of a 
self consistent, non autonomous, delayed differential system: 

p 

= -v°'^{l5 + hr^)+s[^Jc,p{v'^ *h)+I°'*g) (15) 



where 6 is the Dirac function and S is an element- wise function such that S {x") = x°- + + c)aS {x") — b. 



17 



X 



a 



(1 -c)a + h 



{l-c)a + h 



Figure 6. Function S{x°') for the deterministic McKean model given in equation 14 



